Density-matrix renormalization study of the frustrated fermions on the triangular lattice 
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We show that the two-dimensional density-matrix renormahzation analysis is useful to detect the symmetry 
breaking in the fermionic model on a triangular lattice. Under the cylindrical boundary conditions with chemical 
potentials on edge sites, we find that the open edges work as perturbation to select the strongest correlations only 
in the presence of a long range order. We also demonstrate that the ordinary size scaling analysis on the charge 
gap as well as that of the local charge density under this boundary condition could determine the metal-insulator 
phase boundary, which scales almost perfectly with the density of states and the exact solutions in the weak and 
strong coupling region, respectively. 

PACS numbers: 71.10.Hf, 71.27.+a, 71.10.-w 



Geometrically frustrated lattices sometimes provide partic- 
ular situation where the strong correlation between particles 
or spins works effectively to give rise to exotic states. The nu- 
merical approach to such systems made remarkable progress 
in recent years; quantum monte carlo method (QMC) is ap- 
plied to clarify the supersolid state on a triangular and square 
lattices'-^, plaquette state on a quantum dimer model, the ne- 
matic orders in the ring exchanged ladders^. Also the exact 
diagonalization by the symmetry analysis gave much infor- 
mation on the low-energy excitation of the Kagome latticeii, 
nematic orders in the frustrated square lattice^, and so on. In 
the present paper we apply the density matrix renormalization 
group (DMRG) method to the two-dimensional(2D) frustrated 
system. So far, we understand that the systematic treatment to 
2D-DMRG is not known to such system where a competition 
between several orders exists. We give a prototype analysis 
on how to determine the translational symmetry breaking, to 
detect the metal-insulator transition and so on. We choose as a 
representative frustrated system a fermionic system on the tri- 
angular lattice. Since the QMC is quite hopeless to such sys- 
tem due to the minus sign problem, the 2D-DMRG approach 
is useful instead, which can cope with much bigger sizes com- 
pared to the exact diagonalization. However, it is relatively 
quite difficult to get a reasonable results in DMRG when the 
system is critical or in higher dimension than one, which is 
demonstrated in the t-J model on a square lattice^. Recent 
2D-DMRG analysis with the cylindrical boundary condition 
is given by systematically varying the aspect ratio of the finite 
cluster, which successfully determined the magnetization of 
the Neel order in the square lattice Heisenberg model^. How- 
ever, the same analysis for the triangular lattice Heisenberg 
model turned out to be far difficult. In the present paper, we 
adopt the same cylindrical boundary condition, and demon- 
strate that the triangular lattice geometry combined with this 
boundary condition is useful to detect the symmetry breaking 
in an unbiased manner, which is also possible in the case of 
other lattices, e.g., the bipartite square lattice. 

We mention in advance that even though the present analy- 
sis turned out to be successful, there are only few cases where 
the 2D-DMRG is reasonably adopted; (1) in detecting the 
Ising type of long range orders (with quantum fluctuation) 
when the ground state candidates are elucidated, (2) when the 



numerous number of states (more than the order of ~ 10000) 
are kept by e.g., using the parallel computing system^, (3) 
analysis on metallic state away from the phase boundary. Still, 
the difficulty depends on what quantity we measure, and the 
numerically rigorous data provided by the DMRG would be 
of help in many occasions as in the present case which is one 
of the examples of case (1). 

We choose as the simplest fermionic model the t-V model 
whose Hamiltonian given as. 
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Here c| (c^ ) are creation (annihilation) operators of fermions 
and Uj (=c]c^) are number operators. The interactions act 
only between neighboring sites {ij). Anisotropies of the 
hopping amplitudes and repulsion strengths are all positive 
and are given by [tij, Vij) = {t' , V') for the vertical bond 
and (t, V) for the remaining bond directions, as shown in 
Fig.[TJa)-(c). We focus on half-filling where we have compet- 
itive orders due to the commensulability of charges, namely 
one charge per two sites. We take t — 1 as the unit of energy 
and fix i' 1 if not otherwise stated. 

The overall ground state phase diagram of the t-V model is 
presented by part of the authors based on the exact diagonal- 
ization on a 4 X 6 cluster and on a strong coupling analysis^-. 
When the interaction is enough large, the phase diagram is 
characterized by the three different phases according to the 
anisotropy of V and V' as shown schematically in Fig. [Ifd). 
Around the regular triangular geometry of the interactions, 
V' ^ V, we have a partially charge ordered liquid called a 
"pinball liquid". This phase breaks the translational symmetry 
as characterized by the wave vector number, k = ^i^-, 
and it originates from the long range order of "pins", while 
still retains a coherent metallic property. If one enters the 
anisotropic region, V' > V, one finds an insulator where the 
particles align in stripes in the horizontal (or diagonal) direc- 
tion in order to avoid the energy rise by the stronger inter- 
action, V'. We also have similar stripes in V < V which 
extends along the vertical direction. The weakly coupled to 
intermediate coupling region is, however, not clarified yet, 
which will be focused in the present paper. 
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FIG. 1 : (Color online) (a)-(c) Schematic representation of the finite 
cluster, where we take x- and y-direction as OBC and PBC, respec- 
tively. The block states are constructed along the one-dimensional 
chain array shown in arrows. The chemical potentials, /i, are placed 
on each open edge sites. In (a)-(c) the dotted lines which are 
anisotropic bonds with Vij — V' and tij = t' are taken in differ- 
ent directions, (d) Schematic configurations of particles on the clus- 
ter in different phases of the t-V model; (i) horizontal or diagonal 
stripes for V' > V, vertical stripe for V' < V, and pinball liquid for 
V ~ F > 3t. 



The 2D-DMRG calculation is performed on a iV = x 
Ny cluster shown in Fig. [na)-(c). We keep up to 3200 ba- 
sis for each DMRG block and undergo ~ 20 sweeps until the 
ground-state energy converges within an error of ^ lO^^'^t. 
Here, we adopt the open (OBC) and periodic (PBC) bound- 
ary condition in x- and y-direction, respectively as shown 
schematically in Figs, [na)-(c). It is well known that if all 
the boundaries are taken as periodic, the translational symme- 
try of charge density is always preserved at finite system size 
since the wave function is the superposition of all the degener- 
ate states with equal weight. However, if we take the the cylin- 
drical topology of the triangular lattice system and opens part 
of the boundaries, it gives rise to the lifting of the degenerate 
wave functions in the DMRG calculation. Therefore, we can 
easily detect translational-symmetry breaking charge-ordered 
states shown in Fig.fTJd). We additionally place the chemical 
potentials, ^ = {V + V')/2 or V on each open edge site in 
order to suppress the artificial Friedel oscillation induced by 
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FIG. 2: (Color online) Charge density at sites indexed by numbers 
1 to 6 in the cluster Fig. 1(a), where the parameters are taken as (a) 
V = 1, (b) y = V', and (c)V' = 1 under the variation of V or V' . 
The calculations are given on 8 x 6 cluster (Fig. 1(a)) with periodic 
boundary in the y-direction and open boundary in the s-directions. 



OBC. 

Usually, in the exact diagonalization and other finite size 
methods, we figure out what kind of structural factor the state 
has by analyzing the two-point correlation functions. If the 
considered structure has a true long range order, the ampli- 
tude of the structural factor shall remain finite in the thermo- 
dynamic limit after the finite scaling analysis. However, in 
the present 2D-DMRG calculation, each of these phases are 
detected more simply by the spatial structure of the charge 
density as follows; as in Fig.[Tla) we take the PBC on one of 
the l/-bond directions. Then, the other two open bonds with 
different interactions, V and V', artificially lifts the degener- 
acy of the finite size ground state and enables us to distinguish 
several translational symmetry broken and also the unbroken 
phases. Figure |2] shows the charge density of each character- 
istic state for sites indexed by numbers 1 to 6 in Fig.fTla). The 
two-fold periodic structure in Fig.|3a) and (c) indicate the di- 
agonal and vertical stripes, respectively. We also detect the 
three fold periodic rich-rich-poor density of charges (denoted 
as A-A-B structure) in the pinball liquid state in Fig.|2b). In 
contrast to those phases, the metallic phase at small V or V' 
has almost uniform charge density along this periodic bound- 
ary. Therefore, the way we assign the boundary conditions 
should be the one that could efficiently detect the differences 
between ordered states. 

Here, we must note that the horizontal stripe, which is an- 
other possible two-fold periodic structure at V' > V, is not 
realized if we take the boundary conditions as in Fig. \lla). 
Therefore, we consider another type of cluster shown in 
Fig- Ifb) where both the diagonal and horizontal stripes are 
compatible, and compare their energies. The calculations on 
the 6x8 cluster confirm that the horizontal order is always 
the ground state in the insulating phase at atV' > V and the 
diagonal order belongs to a higher state with small excitation 
energy. 

We now present the phase diagram on the plane of V' and 
V in Fig. [3] As mentioned earlier we find three characteris- 
tic phases in the strong interaction region; the pinball liquid 
phase is sandwiched between two different stripes. As has 
been discussed in Ref.|9 10|, the width of the pinball liquid 
phase is determined by the absolute value of t and t', e.g. 
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FIG. 3: (Color online) Phase diagram of the t-V model on the plane 
of V/t and V'/t obtained by the DMRG calculation. The size of the 
circle corresponds to the amplitude of the charge density modulation. 
The calculations are given on 8 x 6 (or 6 x 8) cluster in Fig.[T] 



for t < t'. Since these three phases breaks the translational 
symmetry with different characteristic wave numbers the tran- 
sition between them are of first order. In the strongly interact- 
ing region, the phase boundaries show excellent agreement 
with the one found in the exact diagonalization. The simi- 
lar phase diagram has been obtained by the same 2D-DMRG 
method for t-U -V model' \ and by the variational monte carlo 
(VMC) method'^. The physics of the charge ordering is es- 
sentially the same but the pinball-liquid phase is replaced by 
the three-fold charge-ordered metallic phase with doubly oc- 
cupied sites. 

In one dimension, the open boundary is regarded as impu- 
rities which induces a Friedel oscillation in metals and can 
be used to analyze the wave number and critical exponent of 
the Tomonaga-Luttingerliquid'^ '^. These Friedel oscillations 
are not a true long range order and shall be suppressed if the 
proper chemical potentials are placed on open edges. On the 
other hand, if the system has a true long range order with par- 
ticular symmetry fc ^ 0, the OBC actually lifts the degeneracy 
of wave functions, regardless of whether the chemical poten- 
tials is present or not. We take advantage of this fact in 2D and 
use the partially open edges which would give rise to the spa- 
tial structure of the local quantities (charge density) only when 
we have long range order. The partially opened edge must be 
compatible in geometry with the symmetry of the orders to get 
rid of the the translated counterpart components. We note that 
the chemical potentials not only suppresses the Friedel oscil- 
lation but also helps the variational states not to fall into the 
local minimum in the case of the Ising-type of orders. Even if 
we adopt this method to the square lattice system, such break- 
ing of symmetry will also appear. This is because the transla- 
tional symmetry is broken in a direction with open ends and 
concomitantly the breaking of symmetry in the other direction 
which is periodic occurs. From these results, we may argue 
that the analysis based on the cylindrical boundary condition 
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FIG. 4: (Color online) Finite size scaling of the expectation value of 
(a) charge density extracting the mean charge density at the center 
four sites of the system, and (b) the charge gap at V = V and 
t = t' — 1. The calculations are given on 4 x L cluster with periodic 
boundaries in the i-directions and open boundary in the j/-direction. 



are useful to detect the breaking of translational symmetry in 
the two-dimensional lattices. We also confirmed that the same 
analysis is well adopted to the the Kagome lattice as well'''. 

In the above analysis, we only considered the fixed cluster 
size. However, the phase boundary is usually overestimated 
by the small system size, particularly when we have finite 
but long correlation length compared to the size of the clus- 
ter. Thus one cannot really discriminate the short range order 
from the true long range order without the proper size scaling 
analysis. In order to get the more quantitative information, we 
give the scaling analysis on these local quantities by using the 
L X 4 cluster shown in Fig.[TJc). We focus on the four different 
sites placed at the center of the cluster (see the sites marked 
by numbers 1-4 in Fig.lTfc)) which are the least influenced by 
the open boundary. Figure |4ja) shows the deviation from the 
average charge density of these sites as a function of inverse 
system length in the pinball liquid phase atV = V . Here 
we take L = 6n with integer n, which is a period compatible 
with both the two-fold and three-fold periodicities. As we can 
see, at V — V < 6 the deviation of charge density is extrap- 
olated to zero, while over V' ~ V > 6 they change the slope 
significantly and become finite. This means that the transla- 
tional symmetry is broken at V = V ^ 6, and a three-fold 
periodic long range order appears. This phase boundary is 
smaller than the one estimated from the VMC simulations on 
the same model, Vc ~ 12t^. The VMC estimation based on 
the Fermi sea slater wave function is an upper bound, and the 
present results lie just in between the Hartree-Fock estimation 
and the VMC one. 

We also estimate the charge gap, Ac, of the three-fold pe- 
riodic phase. The charge gap at each system size and 
particle number = N/2 is given as, A{Ne, N^, Ny) = 
E{N, + 1, AT,, Ny) + E[N, - 1, A^„ A^^,) - 2E{N,, N,, Ny), 
where E{Ne, N^, Ny) is the energy of the corresponding clus- 
ter. FigureHfb) shows the extrapolation of A(A^e, A^^;, Ny) to- 
wards Nx = L, L — > cx) with A^j, = 4. We see that the system 
can be regarded as gapless in the thermodynamic limit, which 
is consistent with the previous study'". The estimated phase 
boundary quantitatively agrees with the one by the charge den- 
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FIG. 6: (Color online) (a) The charge gap in the bulk limit at y = 
as a function of V' for the several choices of t'/t. The finite size scal- 
ing is given in the same way as those in Fig.4. (b) t'/t-dependence of 
the phase boundary between the normal weak coupling Fermi liquid 
phase and the horizontal stripe insulator, Vmi- The inset shows the 
t'/i-dependence of t' /Vmi for the same data. 



FIG. 5: (Color online) Finite size scaling of the expectation value 
of (a) charge density after extracting the mean charge density at the 
center four sites of the system, and (b) the charge gap atV = and 
t = t' = 1. In panel (c), the gap and the charge density dispropor- 
tionation in the bulk limit is presented as a function of V' /t. 

sity in Fig. |2] and guarantees the present analysis on the local 
quantities. 

Now we focus on the weak coupling region of the phase di- 
agram which is not well studied yet. In the horizontal stripe 
phase at V' > V, the charge gap opens due to the staggered 
long range order of particles in the anisotropic direction. The 
gap is well estimated by extrapolating the system size in that 
direction. Figure [Sja) shows the ^-dependence of the gap 
for several choices of V' /t with the fixed V = Q. The same 
scaling of the center site charge density is given in Fig.|5lb). 
As plotted in Fig.|5lc), the smooth opening of the extrapolated 
gap suggests the second order transition, while the extrapo- 
lated charge density shows a sharp development at the transi- 
tion point. These are common behaviors in charge ordering 
induced by long-range interactions'^. However, the transition 
point itself shows an excellent agreement between these two 
quantities. We note that the L-dependence of the gap is very 
small in the insulating state, which reflects the localized char- 
acter of the wave function. Resultantly, although the phase 
boundary of Fig. [3] shifts somewhat to a smaller V- or V'- 
region after the scaling, the correction is less than ^ 0.5t. 



Next, by performing the same analysis, we investigate the 
effect of varying the geometry of transfer integrals, t'/t. Fig- 
ure|6ja) shows the extrapolated charge gap for several choices 
of t'/t. By fitting the y'-dependences, the metal-insulator 
(MI) phase boundary is detected. Here, the gap opens linearly 
which makes the estimate precise. The t'/t-dependence of 
V^i is plotted in Fig. |6jb). Just off the non-interacting point, 
V = 0,V' ^ 0, the gap opens since the shape of the Fermi 
surface at half-filling is a regular square and the perfect nest- 
ing takes place. Therefore the obtained insulator is a typical 
charge-density-wave state with a small gap. The value of Vp^j 
significantly changes at < t'/t < 0.5, and then saturates to 
Vp^j ^ 3. At t'/t ^ 1, the system is regarded as the one- 
dimensional chain, where the exact solution of the charge gap 
is given as Vl^^/t' = 2"^. This can be more easily understood 
in the V^^/t-t' /t plot in the inset, where the value of V^-^/t' 
shows a clear crossover from the V"j^j-linear term to the con- 
stant value at around Vp^j 3. 

We now try to understand the characteristic t'/t-linear be- 
havior of V^i/t' in the weak coupling region. Figure [T] shows 
the t'/t-dependence of the non-interacting density of states 
of the anisotropic triangular lattice at the Fermi level at half- 
filling, D{Ef). We find that t' /D{Ef) scales almost linearly 
with t'/t in the wide region, 0.25 < t'/t < 1.5, which is ex- 
actly the region where we find the same linear relation in the 
inset of Fig. |6jb). By interpolating these two characteristic 
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FIG. 7: (Color online) t' / D{Ef) as a function of t' /t at the fixed 
value of t = 1, where D{Ef) is the non-interacting density of 
states of the anisotropic triangular lattice at the Fermi level at half- 
filling. The inset shows the non-interacting density of states for sev- 
eral choices of t' /t. 



relations, we obtain the following relation, 

t' a + (3t'/D{EF) ' 

where a = 0.04, j3 ~ 0.06. The two different linear relations 
found in Fig. |5jc) indicates that the MI transition changes 
its character from the weak to the strong coupling ones, i.e. 
from the charge-density-wave to the charge ordered insulator, 
at t' /t ~ 1.6, which is also the point where the van Hove 
singularity touches the Fermi level at half-filling. 

To summarize we performed the 2D-DMRG analysis on 
the fermionic system on an anisotropic triangular lattice. We 
showed that in the triangular lattice the cylindrical bound- 



ary condition with chemical potentials on open edge site, the 
translational symmetry breaking is detected by measuring the 
local quantity in a reasonable accuracy. Also the finite scaling 
analysis is found to give a further quantitative estimation on 
the phase boundaries. The demonstration is given on the de- 
termination of the metal-insulator phase boundary as a func- 
tion of the degree of the transfer integral, t' /t. It is found 
that at relatively small t' the phase boundary is scaled by the 
non-interacting density of states at the Fermi level, which indi- 
cates that the phase transition is regarded as a charge-density- 
wave transition due to the instability if the Fermi surface. On 
the other hand, at t' > 1.6t, the phase boundary is scaled to 
= 2t', which means that the system undergoes the phase 
transition to the charge ordered state in the one-dimensional 
manner. 

The role of t' /t on the metal-insulator transition has been 
discussed previously in some articles related to organic com- 
pounds based on the mean-field approach'^. Recently, it is 
proved by the uniaxial strain experiment that the MI transi- 
tion temperature of the 6'-ET2CsZn(SCN)4 into the horizontal 
charge ordered state significantly increases when the geome- 
try of the transfer integral varies ast' /t = — 0.52S. The ten- 
dency that the insulator is stabilized due to t' /t is contrary to 
what we get in the present analysis. This might be because the 
present model deals with the half-filled charge density wave, 
while the experiments should be better explained in the 3/4- 
filled extended Hubbard model (which has a different Fermi 
level and a Fermi surface) instead. We expect that even in this 
latter model the Vmi will be scaled by the density of states. 
The present paper shows for the first time that the 2D-DMRG 
provides a reliable way to cope with the difficult parameter re- 
gion where the electronic coiTelation and kinetic energy com- 
pete with each other, and presented the clear crossover from 
the weak to strong coupling. 
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